In this class, we will talk about model explainability but more in the context of data explainability or root cause analysis. In many cases building a very good machine learning model is not an ultimate goal. What is really wanted is the data understanding. A factory wants to know why the product is plagued with a defect, not to predict afterward if there is a defect or not. A football team wants to know which position is the best for scoring a goal, not what's the probability of scoring from a given position. And even when they want a prediction they would love to see the justification to trust the model. Often a nice plot is worth more than sophisticated machine-learning approaches.
import pandas as pd
import numpy as np
import dalex as dx
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
from sklearn.preprocessing import OrdinalEncoder, OneHotEncoder
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.datasets import load_wine
from sklearn.pipeline import Pipeline
from sklearn_extra.cluster import KMedoids
from sklearn.impute import SimpleImputer
from sklearn.model_selection import cross_val_score, StratifiedKFold, GridSearchCV, train_test_split
from sklearn.metrics import recall_score
data = load_wine()
data
df = pd.DataFrame(data['data'], columns=data['feature_names'])
y = data['target']
df['target'] = y
You should already be familiar with many data visualization techniques so we will not train it now. I just want to share a less popular type of data analysis. Usually plotting the target against any feature is not helpful but after some modification, we might be able to see some patterns.
plt.plot(df.flavanoids, y, 'bo')
[<matplotlib.lines.Line2D at 0x1c5c10efcd0>]
For each value, we can plot the average target for data:
Please note that for the line "above that value" the more left we go the higher fraction of data is covered. The same with the "below that value"
for col in df.columns.drop('target'):
tmp = df.sort_values(col)
plt.title(col)
plt.plot(tmp[col], tmp[col].apply(lambda x: tmp[tmp[col] <= x].target.mean()), label="<=")
plt.plot(tmp[col], tmp[col].apply(lambda x: tmp[tmp[col] >= x].target.mean()), label=">=")
plt.plot(tmp[col], np.convolve(np.ones(20)/20, tmp.target, mode='same'), label= "~=~")
plt.legend()
plt.grid()
plt.show()
Ok, let's just train a model. We are not interested in top performance right now so we will skip hyperparameter optimization. Also, we want to find the pattern in the data we have, so we don't split the data into validation and test set.
model = RandomForestRegressor()
x = df.drop('target', axis=1)
y = df.target
model.fit(x, y)
RandomForestRegressor()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
RandomForestRegressor()
plt.plot(df.target, model.predict(x), 'bo')
[<matplotlib.lines.Line2D at 0x1c5c2f4a7d0>]
Dalex is a python package for model explainability. We will use some of its functions to understand the data and the model better. First, we need to create an explainer model. Since we are not interested in checking the model performance but the relation between the data and the target we will use the whole dataset here. In the first case, we might want to use the testing set.
exp = dx.Explainer(model, x, y)
Preparation of a new explainer is initiated -> data : 178 rows 13 cols -> target variable : Parameter 'y' was a pandas.Series. Converted to a numpy.ndarray. -> target variable : 178 values -> model_class : sklearn.ensemble._forest.RandomForestRegressor (default) -> label : Not specified, model's class short name will be used. (default) -> predict function : <function yhat_default at 0x000001C5A6A9D900> will be used (default) -> predict function : Accepts pandas.DataFrame and numpy.ndarray. -> predicted values : min = 0.0, mean = 0.931, max = 2.0 -> model type : regression will be used (default) -> residual function : difference between y and yhat (default) -> residuals : min = -0.33, mean = 0.0073, max = 0.39 -> model_info : package sklearn A new explainer has been created!
fi = exp.model_parts()
The first step will be feature importance. It's a basic analysis where we calculate the global impact of a feature. The idea in dalex default approach is to measure how much the model performance is worsening after removing this feature. Of course, it would require retraining the model, the optimal set of hyperparameters might be different and it might affect the results. To avoid these problems we do not retrain the model. Instead, we simulate its removal by assigning random values to it. To make it more realistic the values are not completely random, we just shuffle this column in a dataframe, do the prediction, check performance and repeat these steps multiple times.
fi.plot()
Another useful tool is a partial dependency plot. For a given feature we observe what's the average output of our model for different values of this feature. For each considered value we set this value for each row in our dataframe and calculate an average prediction.
exp.model_profile().plot()
Calculating ceteris paribus: 100%|█████████████████████████████████████████████████████| 13/13 [00:00<00:00, 20.28it/s]
We can also create similar plots for single rows. Here for each column, we present what would be the output from the model assuming we keep all remaining values and change the value of this one selected feature.
exp.predict_profile(x.iloc[[15]]).plot()
Calculating ceteris paribus: 100%|█████████████████████████████████████████████████████| 13/13 [00:00<00:00, 59.43it/s]
SHAP values are equivalents of Shapley values for the predictive models. It estimates the effect of a particular value of a particular feature for a prediction of a considered row. It's also done by replacing this value with proper sampling and replacing this value and measuring the effect on the prediction.
exp.predict_parts(x.iloc[15], type='shap').plot()
exp.predict_parts(x.iloc[15], type='shap').plot()
The result is based on sampling so the result for the same row can vary
exp.predict_parts(x.iloc[88], type='shap').plot()
exp.predict_parts(x.iloc[88], type='shap').result
| variable | contribution | variable_name | variable_value | sign | label | B | |
|---|---|---|---|---|---|---|---|
| 0 | alcalinity_of_ash = 21.6 | 0.006685 | alcalinity_of_ash | 21.60 | 1.0 | RandomForestRegressor | 1 |
| 1 | alcohol = 11.64 | 0.112978 | alcohol | 11.64 | 1.0 | RandomForestRegressor | 1 |
| 2 | color_intensity = 2.8 | -0.085843 | color_intensity | 2.80 | -1.0 | RandomForestRegressor | 1 |
| 3 | magnesium = 84.0 | -0.002022 | magnesium | 84.00 | -1.0 | RandomForestRegressor | 1 |
| 4 | flavanoids = 1.69 | 0.013539 | flavanoids | 1.69 | 1.0 | RandomForestRegressor | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 8 | malic_acid = 2.06 | 0.003499 | malic_acid | 2.06 | 1.0 | RandomForestRegressor | 0 |
| 9 | proanthocyanins = 1.35 | -0.002735 | proanthocyanins | 1.35 | -1.0 | RandomForestRegressor | 0 |
| 10 | ash = 2.46 | 0.001580 | ash | 2.46 | 1.0 | RandomForestRegressor | 0 |
| 11 | total_phenols = 1.95 | 0.000924 | total_phenols | 1.95 | 1.0 | RandomForestRegressor | 0 |
| 12 | nonflavanoid_phenols = 0.48 | -0.000261 | nonflavanoid_phenols | 0.48 | -1.0 | RandomForestRegressor | 0 |
338 rows × 7 columns
Task For each class find the most representative examples and plot breakdown for them.
def repres_samples(df):
df.dropna(inplace=True)
df_list = [d for _, d in df.groupby(['target'])]
medoids = {}
for d in df_list:
kmedoids = KMedoids(n_clusters=1, random_state=17)
kmedoids.fit(d.drop('target', axis=1))
medoids[d['target'].iloc[0]] = kmedoids.cluster_centers_[0]
return medoids
medoids = repres_samples(df)
for key in medoids:
print(f'Breakdown for the most representative sample of target {key}')
exp.predict_parts(medoids[key], type='shap').plot()
Breakdown for the most representative sample of target 0
Breakdown for the most representative sample of target 1
Breakdown for the most representative sample of target 2
There are other approaches that can be used for model explainability.
Send it to gmiebs@cs.put.poznan.pl within 144 hours after the class is finished. Start the subject of the email with [IR]
Team members:
df = pd.read_csv("DATA.csv")
df.drop(['Patient_ID'], axis=1, inplace=True)
df.rename(columns={'MonkeyPox': 'target'}, inplace=True)
df.head(3)
| Systemic Illness | Rectal Pain | Sore Throat | Penile Oedema | Oral Lesions | Solitary Lesion | Swollen Tonsils | HIV Infection | Sexually Transmitted Infection | target | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | None | False | True | True | True | False | True | False | False | Negative |
| 1 | Fever | True | False | True | True | False | False | True | False | Positive |
| 2 | Fever | False | True | True | False | False | False | True | False | Positive |
df.describe().T
| count | unique | top | freq | |
|---|---|---|---|---|
| Systemic Illness | 24999 | 4 | Fever | 6381 |
| Rectal Pain | 25000 | 2 | False | 12655 |
| Sore Throat | 24998 | 2 | True | 12553 |
| Penile Oedema | 25000 | 2 | True | 12612 |
| Oral Lesions | 24995 | 2 | False | 12512 |
| Solitary Lesion | 24997 | 2 | True | 12525 |
| Swollen Tonsils | 24995 | 2 | True | 12531 |
| HIV Infection | 24999 | 2 | True | 12584 |
| Sexually Transmitted Infection | 24999 | 2 | False | 12553 |
| target | 25000 | 2 | Positive | 15909 |
print(f'Missing values: {sum(df.isna().sum()) + sum(df.isnull().sum())} ')
print(f'Number of duplicates: {df.shape[0] - df.drop_duplicates().shape[0]} ')
Missing values: 36 Number of duplicates: 22944
palette = sns.color_palette('pastel')
fig, ax = plt.subplots(1, 2, figsize=(15, 5))
fig.tight_layout(pad=5)
target_dist = df['target'].value_counts()
ax[0].pie(x=target_dist.values, labels=target_dist.index, colors=palette, autopct='%.0f%%')
ax[0].set_title('Distribution of target class (Monkey Pox)', fontsize=15);
illness_dist = df[df['target'] == 'Positive']['Systemic Illness'].value_counts()
ax[1].pie(x=illness_dist.values, labels=illness_dist.index, colors=palette, autopct='%.0f%%')
ax[1].set_title('Distribution of systemic illness between infected', fontsize=15);
fig, ax = plt.subplots(2, 4, figsize=(20, 8))
ax = ax.flatten()
for idx, feature in enumerate(df.columns.drop(['target', 'Systemic Illness'])):
sns.countplot(x=feature, hue='target', data=df, palette=palette, ax=ax[idx])
ax[idx].set_title(feature, fontsize=20)
ax[idx].set(ylabel=None, xlabel=None)
ax[idx].tick_params(axis='both', labelsize=12)
df_enc = pd.DataFrame(OrdinalEncoder().fit_transform(df), columns=df.columns)
x, y = df_enc.drop('target', axis=1), df_enc['target']
pipe = Pipeline( [ ('Imputer', SimpleImputer(missing_values=np.nan, strategy='most_frequent')), ('RandomForestClassifier', RandomForestClassifier(random_state=17)) ] )
pipe.fit(x, y)
Pipeline(steps=[('Imputer', SimpleImputer(strategy='most_frequent')),
('RandomForestClassifier',
RandomForestClassifier(random_state=17))])In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. Pipeline(steps=[('Imputer', SimpleImputer(strategy='most_frequent')),
('RandomForestClassifier',
RandomForestClassifier(random_state=17))])SimpleImputer(strategy='most_frequent')
RandomForestClassifier(random_state=17)
exp = dx.Explainer(pipe, x, y)
Preparation of a new explainer is initiated -> data : 25000 rows 9 cols -> target variable : Parameter 'y' was a pandas.Series. Converted to a numpy.ndarray. -> target variable : 25000 values -> model_class : sklearn.ensemble._forest.RandomForestClassifier (default) -> label : Not specified, model's class short name will be used. (default) -> predict function : <function yhat_proba_default at 0x000001C5A6A9DB40> will be used (default) -> predict function : Accepts pandas.DataFrame and numpy.ndarray. -> predicted values : min = 0.0, mean = 0.636, max = 1.0 -> model type : classification will be used (default) -> residual function : difference between y and yhat (default) -> residuals : min = -0.971, mean = 4.9e-05, max = 0.963 -> model_info : package sklearn A new explainer has been created!
exp.model_parts().plot()
Systemic Illness has very big impact. Then, Rectal Pain, HIV Infection and Sexually Transmitted Infection, etc.¶exp.model_profile().plot()
Calculating ceteris paribus: 100%|███████████████████████████████████████████████████████| 9/9 [00:02<00:00, 3.58it/s]
Systemic Illness).¶medoids = repres_samples(df_enc)
for key in medoids:
print(f'Breakdown for the most representative sample of target {key}')
exp.predict_parts(medoids[key], type='shap').plot()
Breakdown for the most representative sample of target 0.0
Breakdown for the most representative sample of target 1.0
for key in medoids:
exp.predict_parts(medoids[key], type='break_down_interactions').plot()
The plots above show the contribution of each feature w.r.t. the most representative sample (i.e., completely healthy patient and the infected one).
exp.predict_profile( np.array(list(medoids.values())) ).plot()
Calculating ceteris paribus: 100%|███████████████████████████████████████████████████████| 9/9 [00:00<00:00, 46.17it/s]
score = cross_val_score(pipe, x, y, scoring='roc_auc',
cv=StratifiedKFold(n_splits=5, shuffle=True, random_state=17))
print(f'Roc-auc of the model: {np.mean(score)}')
Roc-auc of the model: 0.6680934276484868
df_enc.dropna(inplace=True)
x, y = df_enc.drop('target', axis=1), df_enc['target']
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=17)
params = {'n_estimators': [10, 50, 100, 150],
'max_features': [3, 6, 9],
'min_samples_leaf': [1, 3, 5, 7],
'max_depth': [3, 5, 10, 15, 20]}
gr = GridSearchCV(estimator=RandomForestClassifier(n_jobs=-1, random_state=17),
param_grid=params,
cv=StratifiedKFold(n_splits=3, shuffle=True, random_state=17),
scoring='recall',
n_jobs=-1)
gr.fit(x_train, y_train)
GridSearchCV(cv=StratifiedKFold(n_splits=3, random_state=17, shuffle=True),
estimator=RandomForestClassifier(n_jobs=-1, random_state=17),
n_jobs=-1,
param_grid={'max_depth': [3, 5, 10, 15, 20],
'max_features': [3, 6, 9],
'min_samples_leaf': [1, 3, 5, 7],
'n_estimators': [10, 50, 100, 150]},
scoring='recall')In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. GridSearchCV(cv=StratifiedKFold(n_splits=3, random_state=17, shuffle=True),
estimator=RandomForestClassifier(n_jobs=-1, random_state=17),
n_jobs=-1,
param_grid={'max_depth': [3, 5, 10, 15, 20],
'max_features': [3, 6, 9],
'min_samples_leaf': [1, 3, 5, 7],
'n_estimators': [10, 50, 100, 150]},
scoring='recall')RandomForestClassifier(n_jobs=-1, random_state=17)
RandomForestClassifier(n_jobs=-1, random_state=17)
gr.best_score_, gr.best_params_
(0.9775011770172717,
{'max_depth': 3,
'max_features': 3,
'min_samples_leaf': 1,
'n_estimators': 100})
recall_score(y_test, gr.best_estimator_.predict(x_test))
0.9768025078369906
maximum depth: 3;
maximum features: 3;
minimum samples leaf: 1;
number of estimators: 100.